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Abstract We studied the problematic of uncertainties 
in the diffuse gamma radiation apparent in stacking anal- 
ysis of EGRET data at low Galactic latitudes. Subse- 
quently, we co-added maps of counts, exposure and dif- 
fuse background, and residuals, in varying numbers for 
different sub-categories of putatively and known source 
populations (like PSRs). Finally we tested for gamma- 
ray excess emission in those maps and attempt to quan- 
tify the systematic biases in such approach. Such kind 
of an analysis will help the classification processes of 
sources and source populations in the GLAST era. 
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, !_i I 1 Introduction 



■ The EGRET era fades away while the Large Area Tele- 
, scope onboard GLAST is in its final stages of hard- 
ware integration. Many unidentified gamma-ray sources 
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are still unidentified, especially at low galactic latitudes, 
where, for instance, not a single supernova remnant could 
be unambigously detected (e.g., [1]). Stacking techniques 
are then a powerful tool to explore, in conjunction with 
spatial cross- localization of sources followed up by Monte 
Carlo analysis (e.g., [2]), if populations as a whole arise 
in the data. The difficulty being, of course, the correct 
handling of the background gamma-ray emission, which 
represent a more difhcult problematic the higher it is. 
This contribution focuses on this very point: devising 
the first ideas for a method of stacking at low Galactic 
latitudes which not only may allow studying already ob- 
tained data (i.e., EGRET), but also be applied to the 
forthcoming observations. 

In the section below we describe the general ideas of 
the stacking technique performed. In section 3 we com- 
ment on the application of this technique to source pop- 
ulation studies in the Galactic Plane (GP), specifically 
in the study of pulsars. Our results are given in section 
4. To verify the implication of them, we performed simu- 
lations, described in section 5. Finally, a short discussion 
is given in section 6. 



2 The stacking technique 

The general stacking method we have applied follows 
that outlined of reference [5]. In order to perform the 
stacking technique and look for a possible collective de- 
tection of gamma-ray emission above 100 MeV near the 
GP, we have extracted rectangular sky maps with the 
selected target objects located at the center. We have 
used EGRET data from April 1991 through September 
1995 —matching the basehne of the Third EGRET Cat- 
alog [4], in galactic coordinates. The extracted maps for 
each particular target were chosen to be 60 x 60 in size, 
in order to have large fields of views and be consistent 
with the EGRET point spread function (PSF). We have 
transformed the coordinates of each map into pseudo- 
coordinates, with the target object at the center. After 
doing this, the maps were co-added, producing the stack- 
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ing. It was also necessary to extract a diffuse background 
map for each target object. For this purpose, we have 
used the diffuse model that is standard in EGRET anal- 
ysis [S]. In order to take into account the existence of 
identified EGRET sources [HE], idealized sources with 
the appropriate fluxes distributed following EGRET's 
PSF as well as the modulated artifacts were added to 
the diffuse background. On the other hand, the uniden- 
tified EGRET sources were not added to the background 
in this study. It was necessary to normalize each one of 
the extracted diffuse maps {Di) for the different expo- 
sures (ci) of the target objects. The extracted diffuse 
map for each target object was also transformed into 
pseudo-coordinates. Finally the diffuse maps for the co- 
added data were obtained as: l/ttotai '^i where Ci are 
the counts diffuse maps (q — eiDi) and etotai — 
To analyze EGRET data we used the standard likeli- 
hood technique based upon gamma ray counts maps that 
were binned in measured gamma-ray energy and spa- 
tially in rectangular projection in Galactic or celestial 
coordinates [7]. The hkelihood function of the EGRET 
data is the probability of the observed EGRET data for a 
specific model of high energy gamma-ray emission, and 
could be written as the product of the probability for 
each pixel: Lg = njPj, where pj is the Poisson probabil- 
ity of observing Uj counts in pixel j when the number 
of counts predicted by the model is 9j . The logarithm of 
the likelihood is used in hypothesis testing and is usually 
more easily calculated. Neglecting the last term (model 
independent) the logarithm of the likelihood is given by: 

logLg^Y.^n,logi9,)-e,]. (1) 
j 

The point-source component of the model consists of an 
"active" source (cq counts located at {aa,Sa)) subject to 
parameter estimation, and "inactive" sources with fixed 
counts at fixed positions. Thus, the total model predic- 
tion for pixel j is given by 

= 9mulGj + 5fciaslO~^i?j+ 

+CaPSF{aa,Sa,j) + Efc P S F {a^ , Sk , j) , ^ ' 

where Ck is the number of counts for the " inactive" source 
at (afe,i5fc); PSF(Q;,5,j) is the fraction of the PSF located 
at {oi,5) that is in pixel j\ Ej is the exposure in pixel j; 
and Gj = EkGkPSF{(j)jk)/EkPSF{<Pjk) (where (0,fe 
is the angle between pixels j and k). The parameters of 
the gamma-ray model are estimated via the Maximum 
likelihood approach. The sum in equation ([1]) is done for 
pixels within an adjustable analysis radius (nominally 
15° for E > 100 MeV). Within this circle, the Galactic 
diffuse radiation model, is scaled by a multiplier, gmuii 
which is estimated by maximum likelihood. Also, a maxi- 
mum likelihood value is used for the level of isotropic dif- 
fuse intensity, gtias ■ This decouples the likelihood point- 
source analysis from uncertainties in the large-scale dif- 
fuse emision model for analyzing point souces in the 
given region of interest. Only the shape of the model 



over the 15° radius circle is used for point-source anal- 
ysis. The expected value of gmui is 1 and gbias is if 
the galactic diffuse model is correct. In order to test the 
significance of a detection, the model of equation @ is 
used in the likelihood ratio test by testing the null hy- 
pothesis, Ca = 0, against the hypothesis that Ca has the 
estimated value where gmui and gbias have their optimal 
values for both hypotheses. This formalism produces a 
"test statistic": TS = —2{lnLQ — lnLi), where Li and Lq 
are likelihood values with and without a possible source. 
TS^^^ is roughly equivalent to the standard deviations. 



3 Stacking technique applied to Galactic Plane 
source population 

The stacking technique for population studies of gamma- 
ray sources has been applied for different class of object 
in several works (for example: radiogalaxies and Seyfert 
galaxies [3], LIRGs and ULIRGs galaxies [5], clusters 
of galaxies [9]). All these studies have been so far ap- 
plied only to high-latitute source populations because of 
the odds to deal with a dominant and structured dif- 
fuse emission. Here, however, we directly step into a new 
methodology to quantify the systematic biases arising 
from significant diffuse contributions. We applied this ap- 
proach in the analysis of EGRET data to study gamma- 
ray emission from pulsars for energies above 100 MeV. 

Pulsars represent astrophysical laboratories for ex- 
treme conditions. Their properties such as densities, tem- 
peratures, velocities, electric potentials, and magnetic 
fields associated with these spinning neutron stars give 
rise to high-energy emission through a variety of mecha- 
nisms. Before the launch of CGRO in 1991 only Crab pul- 
sar (PSR B053H-21), Vela pulsar (PSR B0833-45) and 
Geminga (but not as pulsar in that moment) were known 
as gamma-ray sources. The instruments on CGRO have 
detected a total of 7 pulsars with high significance: Crab, 
B1509-58, Vela, B1706-44, B1951-I-32, Geminga, and B1055- 
52. The weakest (PSR B1951+32) has a statistical prob- 
ability of occurring by chance of 10~^. Not all seven 
are seen at highest energies: PSR B1509-58 is seen only 
up to 10 MeV by COMPTEL and not at 100 MeV by 
EGRET. Sensivity for an individual PSR detection is 
greatly enhanced once studied phase-coherent. The six 
seen by EGRET all show a double peak in their light 
curve. In addition to the six high confidence pulsar de- 
tection above 100 MeV, three additional radio pulsars 
may have been seen by EGRET: B1046-58, B0656-hl4, 
J0218-I-4232. These three all have chance probabilities 
about 5 orders of magnitude less convincing than PSR 
B1951+32. 

More than 1500 radio pulsars are known and it can be 
expected that this number will continue to grow as more 
refined detection equipment is used and spatial coverage 
is expanded. The gamma ray pulsars can be compared to 
other pulsars in term of the derived physical parameters. 
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Fig. 1 Distribution of observed radio pulsars in a period- 
period-derivative diagram derived from ATNF Pulsar Cata- 
logue [11] , The figure is from reference [10] . See text. 

Figure 1 (from reference [lOj) displays the distribution 
of observed radio pulsars in a period-period-derivative 
diagram derived from Australia Telescope National Fa- 
cility (ATNF) Pulsar Catalogue ^llj Q Lines indicating 
the "rotational" age of the pulsars and their dipole field 
strength are also shown as well as the open field line volt- 
age. The gamma-ray pulsars are shown as squares (large 
dark boxes: seven high-confidence gamma-ray pulsars; 
large light boxes: three lower-confidence gamma-ray pul- 
sars). Gamma ray pulsars tend to be concentrated in 
region with high magnetic field (shown by dashed lines) , 
relatively young ages (shown by solid lines) and their 
open field line voltage is high compared to most pulsars 
(dotted lines). Efforts to search for additional pulsars 
in EGRET data have been unsuccessful due to limited 
statistics. In order to apply the stacking technique in the 
study of gamma-ray pulsars we created subclasses of pul- 
sars using ATNF pulsar survey after sorting them with 
different criteria: Surface magnetic flux density [Fb = 
3.2 X 10^^(F/P)^/^, where P is the pulsar's period), best 
estimate of the pulsar distance (13), Spin down energy 
loss rate (E), Energy flux at the Sun (Fe/D^), Spin 
down age (t — P/{2P)). The classes were chosen ac- 
cordingly on what is already known about the gamma- 
ray pulsars that have been detected. We excluded from 
each list of subclasses those pulsars with |6| > 30° and 
the detected EGRET pulsars: Crab, J0633+1746, Vela, 
PSRB1055-52, PSRJ1706-44, PSRB1951-H32 (see for 
example: [ia[l31[10] and also J2229+6114 [Mj). Those 
high-gamma-ray flux pulsars were excluded because they 
inmediately determine the complete stacking problem. 
We first analyzed pulsars individually using the standard 



^ Figure 1 is from 2003, so pulsars discovered since them 
are not included. 



EGRET software. After that, for each class or subclass 
we have generated stacked maps containing N pulsars, 
with N=2, 4, 6, 50. For each stacked map so gener- 
ated, we have then determined the flux, flux error, up- 
per limit, TS, gmui, gbias, gmui/dbias in the center of 
the maps and gmui/gbias averaged over a 6° x 6° box, 
approximately the size of EGRET's PSF for energies > 
100 MeV. 



4 Results 

Figure 2 shows some of the results of our study for all 
the subclasses of pulsars investigated: E, Fb, t, D and 
Fe/D^. The plots represent the TS obtained (left-axis) 
and gmui/gtias averaged over 6° x 6° box (right-axis) 
versus the number of stacked maps. We did not flnd any 
signal in the Fe/D'^ class. Fb and r class both have 
similar behavior with a peak that appears at high N 
values in the sample. On the other hand, E class has a 
peak that is dominated from a few sources. The D class is 
more spread; but a tendency to show high values of TS 
is apparent perhaps for the more nearby sources, then 
fading away when distances of about 1 kpc and larger 
are reached. 

During stacking, we keep continuously track on the 
individual contribution from the diffuse emission model. 
By comparing the individual contribution and evolution 
of the diffuse emission over the growing numbers of sources 
in the stacking sample, we can immediately judge if a 
change in TS is due to a newly added source or rather 
to a odd diffuse emission value/problematic treatment 
during the stacking. There are several cases to distin- 
guish already: (a) steadily accumulating gmui I gbias as in 
Figure 2 d (D) , which points towards an steadly increas- 
ing dominance of the diffuse emission, thus diminishing 
the chance for determine equal TS when adding more 
sources, (b) constant gmui/ gbias as in Figure 2 a, b, c, e - 
which assures that no glitches in the diffuse model have 
an impact on the outcome of the stacking result. Initial 
jumps at the beginning of the stacking are compensated 
in the average after ~ 5-6 sources in sample. 



5 Monte Carlo Simulations 

In order to understand the results obtained above we 
ran Monte Carlo Simulation creating fictitious objects 
with zero gamma-ray flux at random sky positions with 
|6| < 30°, transforming and co-adding the maps, then 
analyzing the stacked map using the same methods as for 
real objects. 1000 simulations were performed for 2, 10, 
20, 30, and 40 objects added with zero gamma-ray flux. 
Examples of the results obtained are shown in Figure 
3, where the TS cumulative distribution is plotted for 
2, 20 and 40 sources added. Accordingly to our results 
there is no more than 3 % chance of obtaining TS > 50 
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Fig. 2 TS{gmui/gbias) versus number of sources added left- 
axis(right-axis) for different subclasses of pulsars investi- 
gated: E, Fb, T, D and Fe/D'^. TS, solid lines; Qmui/gbias, 
dashed lines. See text. 




Fig. 3 Monte Carlo Simulations: Chance probability of TS. 
2 (solid line), 20 (dashed line) and 40 (dot line) sources added 
with zero gamma-ray flux. See text. 



if 2 sources are added; and no more than 6% if 10 (not 
shown in Figure 3), 20, 30 (not shown in Figure 3), or 
40 objects with zero flux are added randomly in the GP. 



6 Discussion 

In this paper a new technique to do source population 
studies using the stacking method in the GP is presented. 
Different subclasses of pulsars were investigated. Among 
all orderings of pulsars, the one sorted by distance (D) 
is the most promising, following the behavior that we 
expected previous to our study: TS signal fading away 
when the diffuse contribution grows. The fact that the 
ordering by Fe/D"^ does not appear to show a signal 
in the stacking is surprising and will be investigated for 
systematics in subsequent studies. The Fb and r classes 
start to show a signal for quite large number of pulsars, 
which may be due to the fact that only individual sources 
affect the stacking . We anticipate that the technique 
explained in this paper will have application in the study 
of Galactic sources (not only for PSRs) in GLAST era. 
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